Astronomy & Astrophysics manuscript no. LsPoisson 


February 2, 2008 


(DOI: will be inserted by hand later) 





o 
o 

(N 



> 

cn 

a^ 
o 

o 
o 



Least-Squares methods with Poissonian noise: an 
analysis and a comparison with the Richardson-Lucy algorithm 

R. Vio\ J. Bardsley^, and W. Wamsteker^ 

^ Chip Computers Consulting s.r.L, Viale Don L. Sturzo 82, S.Liberale di Maroon, 30020 Venice, Italy 

ESA-VILSPA, Apartado 50727, 28080 Madrid, Spain 

e-mail: robertovio@tin.it 
^ Department of Mathematical Sciences, The University of Montana, Missoula, MT 59812-0864, USA. 

e-mail: bardsleyj@mso.umt.edu 
^ ESA-VILSPA, Apartado 50727, 28080 Madrid, Spain 

e-mail: willem.wainsteker@esa.int 

Received ; accepted 



Abstract. It is well-known that the noise associated with the collection of an astronomical image by a CCD camera 
is, in large part, Poissonian. One would expect, therefore, that computational approaches that incorporate this 
a priori information will be more effective than those that do not. The Richardson-Lucy (RL) algorithm, for 
example, can be viewed as a maximum-likelihood (ML) method for image deblurring when the data noise is 
assumed to be Poissonian. Least-squares (LS) approaches, on the other hand, arises from the assumption that 
the noise is Gaussian with fixed variance across pixels, which is rarely accurate. Given this, it is surprising that 
in many cases results obtained using LS techniques are relatively insensitive to whether the noise is Poissonian or 
Gaussian. Furthermore, in the presence of Poisson noise, results obtained using LS techniques are often comparable 
with those obtained by the RL algorithm. We seek an explanation of these phenomena via an examination of the 
regularization properties of particular LS algorithms. In addition, a careful analysis of the RL algorithm yields 
an explanation as to why it is more effective than LS approaches for star-like objects, and why it provides similar 
reconstructions for extended objects. We finish with a convergence analysis of the RL algorithm. Numerical results 
are presented throughout the paper. It is important to stress that the subject treated in this paper is not academic. 
In fact, in comparison with many ML algorithms, the LS algorithms are much easier to use and to implement, 
often provide faster convergence rates, and are much more flexible regarding the incorporation of constraints 
on the solution. Consequently, if little to no improvement is gained in the use of an ML approach over an LS 
algorithm, the latter will often be the preferred approach. 

Key words. Methods: data analysis - Methods: statistical - Techniques: Image processing 

1. Introduction IStarck et al. 1 1200^ . For example, recent wavelet-based 

approache s have been shown to provide excellent results 

The restoration of images is a common problem in iNeelamani et al. 2004b Unfortunately, these 

Astronomy. Astronomical images are blurred due to sev- algorithms are very expensive to implement, prohibiting 

eral factors: the refractive effects of the atmosphere, the ^^^-^ large-scale image restoration problems and on 

diffractive effects of the finite aperture of the telescope, the p.^blems that require the restoration of a large number 

statistical fluctuations inherent m the collection of images -^^^^^^^ Consequently, for many restoration problems, 

by a CCD camera, and instrumental defects. An iUuminat- ^^^^ sophisticated and computationally more efhcient al- 

mg example is represented by the spherical aberrationof gQ^it^ms must be used. In this respect, the algorithms 

th^rimaiy mirror of the Hubble Space Telescope ^WhiteJ ^^^^^ ^ ^^^^^ Least-Squares (LS) methodology repre- 

llMt) that limited the quality of the images before the ^^^^ interesting class. In this paper we discuss two LS 

detector system was refurbished. approaches: direct and iterative. Direct methods, which we 

The widespread interest in this subject has resulted discuss in SecO are the most computationally efficient, 

in the development of a large number of algorithms ^^ile iterative techniques, which we discuss in Sec. E21 

with different degrees of sophistication (for a review, see ^llow for the straightforward incorporation of constraints. 
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In spite of the beneficial characteristics of the LS- 
based algorithms, astronomers typically use techniques 
based on a non-linear approach. Such algorithms are usu- 
ally more difficult to implement, are less flexible and of- 
ten have slow convergence rates. In p articular, the orig- 
inal Richardson-Lucy (RL) algori thm ijRichardson 111973: 
ILucv lll974HShepp k Vardi Il982|) and later modifications 
have seen a great deal of attention. RL can be viewed as 
an Expectation-Maximization (EM) algorithm associated 
with a Poisson statistical noise model. Linear LS meth- 
ods, on the other hand, can be viewed as the maximum- 
likelihood (ML) approach when the noise contaminating 
the image(s) of interest is additive Gaussian with constant 
variance across pixels. For CCD camera noise, the statis- 
tical assumption inherent in the RL algorithm is much 
more accurate than that of the LS approach (see Sect. [21). 
Nonetheless, it is often the case that these two methods 
provide results of similar quality (see Sec.0J. From a par- 
ticular point of view, this fact is disappointing, since it 
means that in certain instances the RL algorithm is not 
able to exploit the a priori statistical information. This is 
particularly relevant when the incorporation of the a priori 
information results in algorithms that are more expensive 
and difficult to implement. 

The aims of the present paper are as follows: I) to de- 
termine the performance of the LS algorithms when the 
noise is predominantly Poissonian; II) to determine when 
the LS and RL approaches will give very similar qualita- 
tive results. We stress that such questions are not only 
of academic interest. In fact, the authors believe that due 
to certain distinct computational advantages, the LS al- 
gorithms should be used whenever their use is warranted. 
On the other hand, we caution that it is not our inten- 
tion to conclude that the LS approach is always the best 
choice. In fact, as we will show, this is not the case. 

In the next section, we present the statistical model for 
CCD camera image formation as well as the approximate 
model that we will use in the paper. After some prelim- 
inary considerations in Sec. 13 in Sec. |2| we will explore 
the convergence properties of two LS approaches. We will 
also discuss the performance of these algorithms on differ- 
ent objects. Finally, in Sec. 01 we will explore in detail the 
convergence properties the RL algorithm. Throughout the 
paper we will present numerical results. Finally, we present 
our conclusion in Section [S] 

2. Statistical Considerations 

Astronomical image data is typically collected with a de- 
vice known as a CCD camera. The fo llowing statistical 
model (see ISnvder et. a,nil99,'^ ligg.'j) applies to image 
data from a CCD detector array: 

b{i,j) = nobj(i,i) + "-bck(j, j) + »T-rdn(j, j)- (1) 

Here, b{i, j) is the data acquired by a readout of the pixels 
of the CCD detector array; i,j = 1,2, ... ,N (without loss 
of generality, square images are considered); nobjihj) is 
the number of object dependent photoelectrons; nbck(*,i) 



is the number of background electrons; and rij-dnihj) 
is the readout noise. The random variables ?T.obj(«i, ji), 
fT-bck(*i, ji), and nrdn(*i, ji) are assumed to be indepen- 
dent of one another and of riobj(i2,i2), "-bck(*2, ^2), and 

"-rdn(j2, ^2) for il, jl 12,. 12- 

For clearness of presentation, it is helpful to use 
matrix-vector notation. We rewrite Eq. as 

b = nobj + ^bck + Wrdn, (2) 

where the vectors have been obtained by column stacking 
the corresponding two-dimensional arrays. The random 
vector riobj has a Poisson distribution with Poisson pa- 
rameter Ax, where x is the true image, or object, and 
A is a matrix that corresponds to the point spread func- 
tion (PSF). Depending upon the assumed boundary con- 
ditions, its structu re is typic ally block-circvilant or block- 
Toeplitz (e.g., see IVio et al. ..200 3: Vogel 2002): n^cfc is 
a Poisson random vector with a fixed positive Poisson pa- 
rameter (3; and Gaussian random vector with 
mean and fixed variance cr^^^-^. 

In the sequel, we will use the following notation to 
denote model 

b = Poisson[Aa;] -f- Poisson[/3 • 1] -I- N(^q ,^2^ ), (3) 

where Poisson[/.i] denotes a Poissonian random vector with 
mean /x, whereas N(^^^cr^) represents a Gaussian ran- 
dom vector with mean /x and variance cr^ (for iid en- 
tries, fj, = ^ and (T^ = CT^). If a^^„ is lar ge, we have 
Poisson [aj^jj^] w iV(^2 ^^2 (see lFeller1ll97ll) . and hence, 
using the independence properties of the random vari- 
ables in Eq. Q we obtain the following approximation 
of Eq. ©: 

b = Poisson[Aa; + 1 • (/3 + a^^J] - af^^. (4) 

In order to simplify the analysis that follows, we as- 
sume the following simplified model 

b = Poisson[Aa;]. (5) 

The analysis is easily extended to the model given by 
Eq. Furthermore, in regions of high intensity Ax 3> 
(3,a'^, in which case model Q is well-approximated by 
model (0. 

A further useful approximation is possible if the ele- 
ments of b are large. In fact, in this case model lO can be 
well approximated with 

b = Ax + z, (6) 

where 2 is a zero mean Gaussian random vector 

2«70 7V(o,i). (7) 

Here symbol " " denotes Hadamard (element- wise) mul- 
tiplication, and 7 = {AxY^^. In other words, through ©, 
the nonlinear noise model jsj can be approximated with a 
linear, additive, nonstationary, Gaussian noise model. Our 
own numerical experiments suggest that bi > 40 is suffi- 
cient. This condition is true in many situations of practical 
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interest (recall that bi is the number of photons detected 
by the ith pixel in the CCD camera). 

From equations © and iQ we see that z has a 
flat spectrum, i.e., the expected modulus of its Discrete 
Fourier Transform (DFT) is constant. However, here, at 
difference with a gaussian white noise process, the various 
Fourier frequencies are not independent from one another 
(e.g., see Fig.^. The reason is that the point- wise multi- 
plication in the spatial domain corresponds to convolution 
in the Fourier domain, and vice versa. Thus, from we 
have 

2(i,j) = [7®iV(o,i)](z,j), (8) 

where the symbol indicates DFT, " (g) " denotes 

convolution, and represents a two-dimensional, dis- 
crete frequency index. Since convolution is a linear opera- 
tion, and {N(^Q i){i, j)} are iid complex Gaussian random 
variables with zero mean and unit variance, {z(z,j)} are 
complex Gaussian random variables with zero mean and 
a constant variance equal to X); j l7(*ij)P- 

3. Performance of the Least-Squares Approach 

The fact that the noise process z has a flat spectrum, pro- 
vides some insight into the performance of the LS deblur- 
ring algorithms in presence of Poissonian noise. In particu- 
lar, since the LS approach corresponds to the assumption 
that the contaminating noise is an additive white noise 
process, the possible worsening of the performance of LS 
algorithms has to be expected due to their inability to take 
into account the dependence between the Fourier frequen- 
cies that characterize the spectrum of z. The question 
then arises, what happens if such dependence is not taken 
into account? The following two arguments suggest that, 
in many astronomical applications, the consequences are 
not so important: 

1. images of astronomical interest have spectra in which 
only the lowest frequencies have intensities that are 
remarkably different from zero. This is a consequence 
of the fact that the PSFs are nearly band limited, 
i.e, they are very smooth functions. Furthermore, if 
a function is in C*^ (i.e., it has k continuous deriva- 
tives) then its spectrum decreases at least as fast as 
l/i/k+i^ Consequently, this constitutes the lower-limit 
with which the spectrum of the images can be expected 
to decrease; 

2. The discrete Pic ard's condition plus the Riemann- 
Lebesgue lemma (jHansen Ill997l) show that the only 
Fourier frequencies useful for the restoration of the im- 
age are (roughly) those where the contribution of the 
object is larger than that of the noise. 

From such considerations, it is possible to conclude that 
in the construction of the deblurred image only a few fre- 
quencies (i.e., the lowest ones) are primarily used, whereas 
most of the remaining frequencies are of only marginal 
importance. For example, in the case of a star-like source 
(i.e., a non-bandlimited function) and Gaussian PSF with 



circular symmetry (a C°° function) and dispersion (in pix- 
els) (Tp, the observed spectrum is again a Gaussian func- 
tion with circular symmetry and dispersion (in pixels) ap 
given by: 

dp « N/{27Tap). (9) 

In several of the numerical experiments presented below, 
we have used N = 256 and ap = 12. In this case, ap « 3.5. 
With the levels used in the simulations, it happens that in 
b{i,j), the noise becomes dominant when (approximately) 
hj > 7j where 10 < 7 < 20 corresponding to different 
noise levels. Hence, the LS algorithms can be expected to 
be almost insensitive to the nature of the noise. 

Another point worth noting is that, although the set of 
frequencies most corrupted by noise are determined both 
by the noise level and by the spectrum of the PSF, it 
is the latter factor that has the most influence. To see 
this we note that, for the case of the star-like source and 
Gaussian PSF considered above, it is possible to show 
that the frequency index (i,j) where the spectrum of 
the signal and that of z have the same level is given by 
r « N^\ii{As/a,)/{7:apV2), where r = (z2 + j2)i/2^ 
is the amplitude of the source, and ct^ is the level of the 
noise spectrum. 

In the next two sections, we will check the reliability of 
the above arguments in the context of two LS algorithms. 

3.1. The Tikhonov Approach 

The linear systems that one encounters when solving im- 
age restoration problems are often highly ill-conditioned. 
Because of this, solutions can be extremely unstable. 
One of the most standard approaches for handling this 
difficulty is known as Tikhonov regularization. In the 
Tikhonov approach, one obtains a stable estimate of the 
solution of the linear system of interest by solving the pe- 
nalized least squares problem 



xx = argmin ( || Ax — b\\l + \'^\\ x II2 



or, equivalently, by solving the linear system 
{A^A + X^I)x = A^b. 



(10) 



(11) 



Here A is a positive scalar known as the regularization 
parameter. The direct solutions of 1)11(1 can be obtained 
very efficiently using fast implementations of the DFT. 
Moreover, there are various reliable and well tested criteria 
that allow for the estimation of A. A standard technique is 
known as the generalized cross-validation (GCV) method. 
With this approach, the optimal value of A is estimated 
via the minimization of the GCV function 



GCV(A) 



\Ax) 



b\\yN' 



[trace{I - A{X))/N^ 



(12) 



Here A. is the matrix that defines the estimator of Ab, 
i.e., Ab — Ax\, and iV^ is the number of pixels in the 
image. For Tikhonov regularization 



AiX) = A{A'^A + A^/)-! A^b. 



(13) 
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It is useful to express model (|10|l and the GCV func- 
tion (|12|l in the Fourier domain: 



|A(z,j)P + A2 



(14) 



and 



GCV(A) ^N^ J2 



|A(*,j)|2+A2 



/ 



iV-1 



E 



(15) 



Here, one can compute both xx and the minimizer of 
GCV(A) very efficiently. 

Figures |21 - El compare the results obtainable with 
this method when noises are stationary Gaussian and 
Poissonian, respectively. The image bp{i,j), contaminated 
by Poissonian noise, has been obtained by simulating a 
nonstationary Poissonian process with local mean given 
by the values of the pixels in the blurred images, i.e. using 
model Four peak signal to noise (S/N) ratios have 
been considered ^: 20, 30, 40 and 60 dB. They corre- 
spond to situations of very low, intermediate, and very 
high noise levels. The PSF used in the simulations is a 
two-dimensional Gaussian function with circular symme- 
try. The image bg{i,j), contaminated by Gaussian noise, 
has been obtained by the addition of a discrete station- 
ary white noise process to the blurred images. Both the 
Gaussian and the Poissonian noises have been simulated 
through a clas sic inverse distribution method (e.g., see 
l,Tohnsonlll987i) by transforming the same set of random 
uniform deviates. They have exactly the same variance. 
Here, the subject of interest is superimposed to a sky 
whose intensity, in the blurred image, is set to 1% of the 
maximum value of the image. This means that, contrary to 
the Gaussian case where the noise level is constant across 
the image, in the Poissonian case the noise is mostly con- 
centrated in the pixels with highest values. In spite of 
this fact, these figures show that the results provided by 
Tikhonov coupled with GCV are quite similar regardless 
of whether the noise is of Gaussian of Poissonian type. 

These results can be explained if one considers 
Eq. (|14|l , where it is clear that the role of A is to replace the 
Fourier coefficients A{i,j) with small modulus, i.e., those 
coefficients that make the deblurring operation unstable. 
According to the two points mentioned above, the "opti- 
mal" value of A should replace all the Fourier coefficients 
whose modulus is smaller than the expected level of the 
noise in the Fourier domain. Since in bp{i,j) and bg{i,j) 
the level of the noise is the same, such a replacement will 
be quite important and will have a similar effect for both 
images. This is shown by the results in Figs.|Hl-^2 In Par- 
ticular, the c) panels show that GCV chooses A so that the 

^ Here S/N = 20 log(&,nax/&mfx) dB, where bmax is the max- 
imum value in the image Ax. 



frequencies corresponding to the flat part of the spectra 
(i.e., those dominated by the noise) are filtered out. The 
consequence of this is that, for both Gaussian and the 
Poissonian noises, almost the same number of coefficients 
are filtered. Moreover, as is shown in the d) panels, the 
coefficients |&p(«, j)| and |6g(i,j)| corresponding to the co- 
efficients A(i, j) with the largest modulus, are very similar. 
From this, one can conclude that the deblurred images x\ 
can be expected to be very similar regardless of the nature 
of the noise. 

The reason why the two GCV curves are almost iden- 
tical (see the b) panels) is that, independently from the 
nature of the noise, in Eq. (|15() the quantity b{i,j) can be 

replaced by b{ij) = A{iJ)x{i,j) + Zn{i,j), where Zn{i,j) 
is given by Eq. ||SJ| or by a stationary white noise process. 
Now, taking the expected value of the resulting GCV(A), 
it is not difficult to show that 



E[GCV(A)] ^N^ 



\A{i,mh])? 



i.J=0 



(|A(Z,J)|2 + A2)2 



/ 



^JA(*,j)P + A2 



(16) 



Since the term cr^^ is constant, the E[GCV(A)] function is 
independent of the specific nature of the noise. The same 
is not true for the variance. However, because of the argu- 
ments presented above, no instabilities are to be expected. 
This is supported by the fact that in our numerical exper- 
iments we have never experienced stability problems (see 
also Fig.El). 

3.2. The iterative approach to regularization 

Iterative algorithms are commonly used in deblurring 
problems. Although, computationally less efficient than 
the direct methods, such as the Tikhonov approach dis- 
cussed above, they are much more flexible in that they 
allow for the straightforward incorporation of constraints. 
These algorithms provide regularization via a semiconver- 
gence property; that is, the iterates first reconstruct the 
low frequency components of the signal, i.e. those less con- 
taminated by noise, and then the high frequency ones. In 
other words, the number of iterations plays the same role 
as the regularization parameter A in the Tikhonov ap- 
proach. 

Semiconvergence has been rigorously proved only for a 
limited number of algorithms. For others, some theoretical 
results are available but, the primary evidence stems from 
many years of success in use on applied problems. 

The prototypical iterative algorithm for least squares 
problems is the Landweber method (LW). If is the start- 
ing image (often — Q), then the iterations take the form 

Xk = Xk-i + ujA^ [b - Axk-i] , (17) 

where, k — 1,2,..., and cj is a real positive parameter 
satisfying < w < 2/|| The values of u) determine, 
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in part, the convergence of the iteration. The semicon- 
vergence property of this algorithm is typically proved 
using arguments based on the singular values de compo- 
sition of the matrix A (for a discussion of this, see lVogefl 
120021) . However, it is, perhaps, more instructive to rewrite 
Eq. (|17|l in the Fourier domain, obtaining 



2\k 



l-{l-u:\A{^,3r) 



(18) 



with < < 2/ max |^(i, j)| 



If, as usual, the PSF is 



= 1 



assumed to have unit volume, then max 
and < a; < 2. From this equation, one can see that, 
for a given frequency index the closer the term 

a;|j4(z,j)P is to one, the more rapid is the convergence 
to b{i, j)/ A{i, j), which corresponds to the unregularized 
solution. Since, as mentioned above, the largest values of 
the spectrum of A{i,j) correspond to the lowest frequen- 
cies, it is evident from (|18|l that the lower frequencies are 
restored in early iterations, while progressively higher fre- 
quencies are restored as the iteration progresses. 

3.2.1. Convergence properties 

Equation (|17|l shows that the convergence of LW is driven 
by the rate with which the term 



2\k 



ICkii,j) = {l-u;\A{z,j)\') 



(19) 



goes to zero. In order to understand what this means in 
practical situations, it is useful to see what happens in 
the case of a noise-free image when the PSF is a two- 
dimensional Gaussian with circular symmetry and disper- 
sion (Tp. Without loss of generality, we set oj = 1. Then 

£fe(*,j)« [l-exp(-rV?^)]\ (20) 

where = + and ap is the dispersion of the PSF in 
the frequency domain (see Eq. lO). From this equation, it 
is not difficult to see that, even in case of moderate values 
of k, the term within square brackets on the rhs of Eq. (|18|l 
can be well-approximated by the Boxcar function 



1 if 0< |r| <ro.5,fc; 
otherwise 



(21) 



where ro.5,/c is the value of r for which ICk {r) =0.5 (see also 
Fig. [TJ. Therefore, the iterate (|18|l can be approximated 

by 



Xk(l,j) = -^——IlkihJ)- 



(22) 



The requirement that K,k{i,j) < e, with < e < 1 
implies 

In [1 - exp(-r"^/cr^)J 

This result shows that the restoration of the highest fre- 
quencies (i.e., large r), requires a number of iterations that 



becomes rapidly huge. For the case of the star-like sources, 
where all the frequencies have to be restored, this means a 
terribly slow convergence. More specifically, from Eq. H23|l 
one can see that for frequencies (?, j) such as r ^ ct^. A; (x r, 
while for larger values of r, /c increases exponentially. For 
example, some of the experiments presented in this paper 
are based on images with size 256 x 256 pixels and with 
a Gaussian PSF with ap w 3.5. In this case, if e = 0.5, 
in order to have ro.5,fc = ap,2ap,'Sap,Qdp (i.e., ro.5,fc = 
3.5, 7, 10.5, 21), it is necessary that A: « 2, 4, 5600, 3 x 10^^, 
respectively. 

The obvious conclusion is that LW is useful only for 
the restoration of objects for which the low-frequencies are 
dominant, e.g. extended objects with smooth light distri- 
butions. 



3.2.2. Numerical results 

Since, regardless of the noise type, LW reconstructs the 
lower frequency components of the image first (i.e., the 
frequencies where the contribution of the noise is negligi- 
ble), we expect the following for both bp{i,j) and bg{i,j): 

1. the resulting deblurred images should be very similar; 

2. in early iterations the convergence rate of the algo- 
rithms should be almost identical. 



These statements are supported in Figs. I12I15I and 
Fiefs. llfilTi^ In particular, from the last set of figures one 
can see that the convergence curves are almost identical 
until the minimum rms of the true residual is reached. 
After that, because the high frequencies (the ones that 
are more sensitive to the nature of the noise) begin to be 
included in the restoration, the curves diverge. 

4. Richardson-Lucy Algorithm 

In the previous sections we have shown that the LS meth- 
ods are relatively insensitive to the specific nature of the 
noise. However, this does not mean that they are optimal. 
In principle, methods that exploit the a priori knowledge 
of the statistical characteristic of the noise should be able 
to provide superior results. 

In particular, model jSJ) motivates the use of the 
Richardson-Lucy (RL) algorithm for estimating x. RL can 
be viewed as the EM algorithm corresponding to the sta- 
tistical noise model (jSJ. The RL algorithm is defined by 
the iteration 

Xk+i=XkQA^ — , (24) 

Ok 

where bk — Axk, and the fraction of two vectors denote 
Hadamard (component-wise) division. 

Since RL exploits the a priori knowledge regarding the 
statistics of photon counts, it should be expected to yield 
more accurate reconstructions than an approach that does 
not use this information. In reality, as shown by Figs. \2Ul 
1771 the situation is not so clear. These figures provide the 
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convergence rates and the performances of RL and LW 
methods for objects with a size that is increasing with re- 
spect to the size of the PSF (a two-dimensional Gaussian 
with circular symmetry). Two different types of objects 
are considered: a two-dimensional Gaussian and a rect- 
angular function. Since the first target object presents an 
almost band-limited spectrum, whereas for the second tar- 
get object the high-frequency Fourier components are im- 
portant, their restorations represent very different prob- 
lems. For both experiments, a background with an inten- 
sity of 1% of the maximum value in the blurred image 
has been added. Finally, two different levels of noise have 
been considered corresponding to a peak S/N of 30 and 
40 dB, respectively. The first case provides a background 
with an expected number of counts approximately equal 
to 30, i.e., a level for which the Gaussian approximation 
of the Poissonian distribution is not very good. 

From Figs. I20I27I it appears that the performance of 
RL for objects narrower than the PSF is, in general supe- 
rior to LW for the band- limited target. The same is not 
true for the other objects. Interestingly, though, for ex- 
tended objects, i.e. smooth objects with high intensity pro- 
files over large regions, the performance of RL is roughly 
equivalent to that of LS (to properly compare the conver- 
gence rate, it is necessary to keep into account that, for 
each iteration, RL requires the computation of twice the 
number of two-dimensional DFT than is required by LW) . 
This is especially true for the images characterized by the 
best S/N. Motivated by these numerical results, we seek 
answers to the following questions: (i) why does RL per- 
form better than LS on star-like objects, and (ii) why do 
the RL and LS approaches yield roughly the same results 
on extended objects? 

4.1. RL vs LS: preliminary comments 

It is important to note that in practice, when either the RL 
or LS approaches are used in solving image reconstruction 
problems, the exact computation of the maximum likeli- 
hood estimate (MLE) is not sought. For example, as was 
stated above, the LW iteration implements regularization 
via the iteration count. In fact, the objective when using 
LW is to stop the iteration late enough so that an accurate 
reconstruction is obtained, but before the reconstruction is 
corrupted by the noise in the high frequency components 
of the image. Notice, for example, that in Figs. I2DE3 the 
relative error begins to increase at a certain point in both 
the RL and LW iterations. 

As was stated above, one can show that the LW it- 
erates are regularized solutions of Ax = A^b via 
the singular value decomposition (SVD) of the matrix A. 
Unfortunately, such an analysis of RL is impossible due 
to the nonlinearity in the RL algorithm. In particular, 
note the Hadamard multiplication and division in algo- 
rithm (|17l) . Instead, we first note that if A is an invertible 
matrix, RL iterates converge to t he MLE c orresponding 
to the statistical model ((SJ (see IWu Ill98.^ . Hence, RL 



can be viewed as an EM algorithm l)Carasso 1119991) . The 
MLE is also the minimizer of the negative log-likelihood 

function associated with Eq. ij^l; namely, 

J{x) = l^[Ax-bQlog{Ax)]. (25) 

When the RL algorithm is used on image deblurring 
problems it exhibits a similar convergence behavior to 
that of the LW iteration. Specifically, for ill-conditioned 
problems, the RL iterates {xk} provide more accurate re- 
constructions in early iterations (semiconvergcnce p rop- 
erty) , while in later iterations blow-up occurs ( Lucv 1974i 
ICarasso1ll999|) . To explain why this occurs, we note that 
the function J is convex. In fact, assuming A is positive 
definite, J is strictly convex. In this case, the minimum 
of J is the unique solution of the equation VJ(a;) = 0, 
where V J is the gradient J; that is. 

It is clear that provided A is invertible, the solution of 
(|26|l is obtained when Ax = b. That is, x^ = A^b is 
the solution of Eq. H26|l . Under the same conditions on 
A, the LW iteration converges to the same value. Thus, 
we would expect that since blow-up occurs as Xk — s- a;* 
in the LW iteration when A is a poorly conditioned ma- 
trix, we will see the same results when we use RL. One 
consequence of this fact is that during reconstruction, RL 
uses, effectively, only a few frequencies and therefore can- 
not fully exploit prior statistical information regarding the 
noise (this should require the use of the entire spectrum). 

4.2. RL vs. LS: a sensitivity analysis 

In order to obtain a deeper understanding of the RL and 
LS approaches and to answer the two questions posed 
above, it is useful to introduce the quantities 

ALw(a;fc) - -A^rk-, (27) 
AnL{xk)^-XkQA^^, (28) 

which provide the correction to the solution Xk at the k-th 
iteration for LW and RL, respectively. Here, 

rfc = Axk - b, (29) 

and, without loss of generality, we have set w = 1 in the 
LW iteration. We note that in order to obtain Eq. H28|l 
we needed the identity A^l = 1 to hold. From these 
equations it is evident that at each iteration LW corrects 
the solution Xk with a quantity proportional to r^, while 
the correction provided by RL is proportional to x^ it- 
self. Thus it is not surprising that RL outperforms LW 
when applied to reconstructing objects composed of star- 
like sources on a flat back ground, since in the early stages 
of both iterations the entries of Xk are large and increas- 
ing in regions corresponding to the positions of the objects 
in the image, while the values of are correspondingly 
small and decreasing. 
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However, as has been shown by the simulations pre- 
sented above, RL does not outperform LW when applied 
to reconstructing objects with smooth light distribution 
and whose spatial extension is broader than the PSF. In 
order to understand this phenomena, it is useful consider 
the negative Jacobian matrices of the quantities H27|) and 
(EHI): 

JLw{xk) = A^A (30) 



and 



jRL{xk) = diag 



Axk 

T 



+ diag[a;/j] A diag 



Axk Axk 



(31) 



These matrices provide the sensitivities of the LW and RL 
algorithms to changes in the components of the iterate x^- 
Equations H30|l and (|31() allow us to make several obser- 
vations. We begin by considering (|3DJ). Since in general 
astronomical applications A is the discretization of a PSF 
with an almost limited spatial support, the LW sensitivity 
matrix A^ A will also have spatial support that is almost 
limited. From this observation, we can conclude that for 
a given pixel the corresponding component of the vector 
Alw will be most sensitive to changes in the value of the 
pixel itself and in the values of "nearby" pixels. Here, the 
term "nearby" is defined by the characteristics of the PSF. 
More specifically, as the spread of the PSF increases, so 
does the collection of "nearby" pixels. 

Perhaps an even more important observation, is that 
the sensitivity of the LW iteration to perturbations in a;^ 
is independent of both Xk and b. Consequently, the al- 
gorithm has no means of distinguishing between low and 
high intensity regions within the object, and hence, per- 
turbations of the same magnitude are allowed for compo- 
nents of Xk corresponding to regions of both low and high 
light intensity. This explains why in areas of low light in- 
tensity (where the margin of error is very small) LW, and 
the least squares approach in general, does poorly. 

The sensitivity matrix H31() for the RL iteration is more 
difficult to analyze. However, some simplification is possi- 
ble when one considers the problem of the restoration of 
a flat background or of regions of an image in which the 
intensity distribution varies smoothly (e.g., the interior of 
the rectangular function considered in the simulations). In 
fact, in this case it is possible to define a region f2 where 
the image can be considered constant or almost constant. 
Because of the semiconvergence property of RL, in such 
regions the components of the vector converge rapidly 
to zero (this has been verified via numerical simulations). 
Thus the first term in l|31|l converges to zero. The same 
is not true for the second term. Thus, it is reasonable to 
expect that it will provide an accurate approximation of 
the sensitivity of RL within fJ. 

Provided that the spread of the PSF is small relative 
to the size of fi, early in RL iterations the vector x^ is ap- 
proximately constant and close to b, i.e. those pixels val- 
ues are reconstructed rapidly, within Q.. Hence, the vector 



b/{Axk Axk) ~ 1/ Axk is also approximately constant 
within Q.. In addition, if we define Dq, to be the diagonal 
matrix with components 



\D 



I. j en 
0, Jin 



then 



DnAxk ~ DnADnXk 



(32) 



(33) 



will be accurate within the interior of fi. To obtain H33|l 
we used the fact that Xk is approximately constant on n 
and that the spread of A is small compared to the size of 
n. Finally, the second term of l|31|l can be approximated 
within n as follows: 

Dn<^\a.g[xk\A^dia.g[b/{Axk Axk)\A 

« di&g[xk\DnA^ DnAi&g[b/ {Axk Axk)]A (34) 
« di&g[Dn{xk/Axk)]A^A (35) 



DnA'A. 



(36) 



Approximation (|34|l follows from H33|) . Approximation 
(j35|l follows from the fact that, as stated above, early in 
RL iterations b/{Axk Axk) « 1/ Axk is approximately 
constant. Thus we see that not only does the second term 
in 131|) not converge to zero, it is well-approximated within 
n by the LW sensitivity (|30(l . Recalling that the first term 
in H31|l converges rapidly to zero in RL iterations, it is 
therefore not surprising that RL and LW provide similar 
results in the interior of the rectangular object mentioned 
above. We can extend this discussion to extended objects 
in general by noting that such objects can be decomposed 
into a union of regions in which the light intensity is ap- 
proximately constant. Hence, RL and LW should provide 
similar results for extended objects in general. 

4.3. RL vs. LS: convergence properties 

As shown above, LW presents an acceptable convergence 
rate only in case of restoration of extended objects. 
Unfortunately, understanding the convergence properties 
of the RL algorithm (|17() is a very difficult affair since it is 
not possible to carry out an analysis similar to that done in 
Sect. 13.21 For this reason, in order to obtain some insight, 
we consider, again, a noise-free signal b and Gaussian PSF 
with circular symmetry and variance ct^. In addition, we 
suppose that the object of interest is a circular Gaussian 
source with variance cr^. The amplitude of the source is 
not considered since RL conserves the total number of 
counts. Due to the connection between the RL and LW it- 
erations discussed in the previous section, an understand- 
ing of RL convergence may provide further understanding 
of the convergence of the LW iteration. For simplicity, in 
what follows we work in the continuous, and results will 
be later discretized. 
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If a Gaussian function with variance is denoted by 
G'[(T^], the foUowing facts are useful: 

G[al] 



G 



9 9 



G[al] G[al] = G 



G[<jI] (g) G[cr|] ^G[al + al 
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(37) 

(38) 
(39) 

Here, the symbol " 0" indicates convolution. From these 
equations it is evident that the result of any of the above 
operations produces a new Gaussian function. Only the 
first operation requires a condition be satisfied, i.e., cr| > 
a\. This will always be satisfied during the RL iteration 
(see Eq. (1221) below). 

A point worth noting is that if we define to be the 
variance of the Gaussians on the right hand side of (|37|) - 
H39(l . then for equation H37() we have a\ < cr^; in equation 
1)38(1 we have < a\ < a\\ and in equation H39|l we have 
(T^ — af+(T2 - Consequently, only the operation l(SSI) results 
in a Gaussian function with a variance that is smaller than 
both ai and a2- 

If the true object a; is a Gaussian with variance a^, 
then using (|39(l and the fact that Ax = b, it is 



^2 2 , 2 



(40) 

As an initial guess in the RL algorithm, we take Xq — 1. 

Now, let's suppose that Xk = G[cr^]. Using Eqs. (pTzjl 
- ()39(l one can obtain 



G 



It is a straightforward exercise to show that 



k+l 



al + alal 



2\2 



< 1 



(41) 



(42) 



provided 

al>al-al. (43) 

To prove that (gSl holds for all k, we use induction. First, 
note that since Xq = 1, A^l = 1, and A'^ = A, we 
have that Xi = A^b = G[ap + a^]. Then aj = + a^, 
and hence, Eq. (|43|l is satisfied for k = 1. Now, we show 
that if Eq. 14311 holds for k, it must hold also for k + l. 
By replacing (t^_|^]^ in Eq. (|43|l by the argument of the 
Gaussian function on the right hand side of Eq. 141f) . one 
can obtain an equivalent inequality involving cr^ given by 

q{al) > 0, (44) 

where q{cr'^) defined by 

qia') = a'ia'al + ^ + a'a^ + (a^ - al)ia' + a^f 
= 2al{aY + {ial 2alal)a^ + (a^ - ala% 

Notice that q is a quadratic function. We can therefore 
find its zeros via the quadratic formula. These are given 

by 



' ^2 2 



(45) 



Since ap > 0, we know that the graph of q is an up- 
ward opening parabola. Furthermore, by H40|l we have 
cr^ — (7p = > 0, and hence, we know that if > af — ap, 
then g(fT^) > 0. Thus 1)44(1 follows from the inductive hy- 
pothesis, and our proof is complete. 

In light of these findings, it is possible to consider some 
convergence properties of the RL algorithm. We begin by 
showing that the sequence {cr^} converges to — a^ — ap. 
First, note that Eqs. and (|^ imply that {cr^} is a de- 
creasing sequence that is bounded below by a^—ffp. Hence, 
{cr^} converges to some a* > al — ap. From inequalities 
(|^ and 1031) I we have that R{al) = 1. Furthermore, 
the arguments used in the proof of ((43(1 imply that if 
CT^ > cr^ — (7p then R{(T^) < 1. Thus it must be that 

9 2 2 

In regard to the convergence rate of the RL algorithm, 
Eq. 1(42(1 shows that, almost independently from the char- 
acterists of the object, in the very first iteration, when 
(Tf, ^ (Tp, we have 



^2/^2 



0. 



(46) 



In fact, if a;o = 1 (i.e., CTq = 00), it is not difficult to 
see that Xi — Ab, i.e., the result of the first iteration is 
given by G{ap + u^). At this point, there are two possible 
situations: 



CT? > CTp. In this 



1. For extended objects, we have cr^ ~ "1 

case, R{<j\) ~ 1. In general, this means that we can 
expect rapid progress in early iterations; after that the 
convergence rate slows down remarkably. This behav- 
ior is similar to that of the LW algorithm; 



2. For stellar-like objects, we have a 
set Ok — a Up, then 



(l + a2)2 



ai. Now, if we 



(47) 



For example, when a = 1, 0.5, 1/3, 1/6, then i?(cr^) = 
0.750, 0.960, 0.990, 0.999, respectively. In other words, 
although the convergence rate of RL slows down as the 
iteration progresses, this effect is not as pronounced as 
it is for the LW algorithm. 

These statements are confirmed by Figs. |2H1 1221 

A comparison between the RL solution for star-like 
sources at the k-th iterate 



Xfc(i, j) cx exp(-r2/2CT^) 
with the corresponding LW solution 
Xfe(i,j) = nfe(i, j), 



(48) 



(49) 



provides some additional insight into the convergence 
properties of these algorithms. In fact, from Eq. (I48II it 
is evident that although the high frequencies are filtered 
in the RL algorithm, the filter is less stringent for high 
frequencies than is the Landweber filter. The consequence 
is that, in general, at a given fc, RL has available a broader 
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range of frequencies to restore the object. On the one 
hand, this can improve the convergence rate of RL com- 
pared to LW; on the other hand this could create problems 
when one or more star-like objects are superimposed with 
an extended object. In fact, a few RL iterations are suf- 
ficient to restore the extended object. The same is not 
true for the star-like objects. Therefore, more iterations 
are necessary. However, because of the amplification of 
the noise, this means the degradation of the results in the 
parts of the image not occupied by the star-like objects. 

This effect is clearly visible in the experiment shown 
in Fig. EOI 

5. Conclusions 

In this paper we provide explanations for why, in spite 
the incorporation of a priori information regarding the 
noise statistics of image formation, the RL deblurring al- 
gorithm often does not provide results that are superior 
to those obtained by techniques based on an LS approach. 
In particular, we have identified a possible explanation in 
the regularization approaches of the specific algorithms. 
In fact, the adoption of a priori smoothness constraint 
in the Tikhonov approach, or the need to stop the iter- 
ations before blow-up occurs in the iterative approaches, 
e.g. both LW and RL, do not permit the full exploitation 
of the information contained in the highest Fourier fre- 
quencies, i.e., those where the specific nature of the noise 
has the largest influence. This has two consequences: I) 
the performance of the LS algorithms is almost insensi- 
tive to whether noise is Gaussian or Poissonian; II) the 
RL algorithm does not fully benefit from the fact that it 
incorporates the specific statistical model of the noise. In 
other words, the regularization of the solution implies a 
levelling out of the possible performances. In this respect, 
much more than a detailed knowledge of the nature of 
the noise is needed. Specifically, some rough a priori in- 
formation regarding the solution, e.g. is it an extended or 
star-like object, is needed before one can know whether or 
not RL will provide superior results. Our numerical exper- 
iments support these conclusions. In particular, the fact 
that reconstructions gotten via the RL algorithm are often 
comparable to those of LW, i.e., an unsophisticated and 
very slow algorithm, indicates that resorting to advanced 
and often complex techniques is not always justified. 

We stress that such conclusions are not only of aca- 
demic interest. In fact, with respect to the ML algorithms, 
in general the LS algorithms are much easier to imple- 
ment, are more flexible as concerns the incorporation of 
constraints, are more amenable to a theoretical analysis of 
their characteristics and are computationally less costly. 
Consequently, unless the use of a different approach is jus- 
tified, they should be considered the standard approach. 
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Fig. 1. Comparison of Poissonian vs. stationary Gaussian 
noise and corresponding power-spectra for the satellite im- 
age shown in Fig. 0] The variance of the two noises is the 
same. 
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FoEson /rm^= 13 - rmsj = 33 





Foisaon/rrm^= 121 - ■rrtis; = 322 





Blurred lm^« / i3a 




Gau^an / rriE ^ = 13 - rm&, = 32 



Gau&&an/ rrn& ^= 122 - rm&,= 320 





Fig. 2. Comparison of the results obtained by Tikhonov 
coupled with GCV in case of Poissonian and Gaussian 

noises (see text). The images have sizes 256 x 256 pixels, 
the PSF is Gaussian with circular symmetry and disper- 
sion set to 12 pixels, S/N = 20 dB. rmsT and rmsj are 
the rms of the true residuals calculated on the entire im- 
age and only on the pixels corresponding to the satellite, 
respectively. 



Fig. 3. Comparison of the results obtained by Tikhonov 
coupled with GCV in case of Poissonian and Gaussian 

noises (see text). The images have sizes 256 x 256 pixels, 
the PSF is Gaussian with circular symmetry and disper- 
sion set to 12 pixels, S/N = 30 dB. rmsy and rmsj are 
the rms of the true residuals calculated on the entire im- 
age and only on the pixels corresponding to the satellite, 
respectively. 
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Fig. 4. Comparison of the results obtained by Tikhonov 
coupled with GCV in case of Poissonian and Gaussian 

noises (see text). The images have sizes 256 x 256 pixels, 
the PSF is Gaussian with circular symmetry and disper- 
sion set to 12 pixels, S/N = 40 dB. rmsT and rmsj are 
the rms of the true residuals calculated on the entire im- 
age and only on the pixels corresponding to the satellite, 
respectively. 



Fig. 5. Comparison of the results obtained by Tikhonov 
coupled with GCV in case of Poissonian and Gaussian 

noises (see text). The images have sizes 256 x 256 pixels, 
the PSF is Gaussian with circular symmetry and disper- 
sion set to 12 pixels, S/N = 60 dB. rmsy and rmsj are 
the rms of the true residuals calculated on the entire im- 
age and only on the pixels corresponding to the satellite, 
respectively. 
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Fig. 6. Relationship between the estimates Xg and Xp, 
corresponding to the Gaussian and the Poissonian noise, 
respectively, obtained in 100 different realizations of the 
experiment shown in Fig.|3| The mean value and the stan- 
dard deviation are 0.128 and 4.55 x 10^* for Xg, and 0.128 
and 9.10 X lO""^ for A„. 




Fig. 7. Plot of 1 — JCk{r) vs. r for different values of the 
iteration count k. 



14 



R. Vio, J. Bardsley, & W. Wamsteker; Algorithms for Image Restoration 



Panel a) 




Panel b) 



6.86 
> 6.84 
6.82 




0.005 0.01 0.015 0.02 0.025 0.03 



Panel d) 



0.05 

0.04 

S 0.03 
< 

0.02 
0.01 



Fig. 8. Comparison of the results obtained by Tikhonov 
coupled with GCV in case of Poissonian and Gaussian 
noises, This figures correspond to the experiment shown 
in Fig. 12 Panel a) the coefHcients \bg{i,j)\ and \bp{i,j)\ 
in decreasing order. The two horizontal lines represent the 
values of A for the two noises; b) corresponding GCV func- 
tions; c) coefficients \bp{i,j)\ and \bg{i,j)\ corresponding 
to the first 2000 coefficients of j)| shown in panel a). 
The vertical lines show the indices of |&(i,j)| closest to 
A. d) Ab{i,j) = \bg{i,j) - bp{i, j)\/\bg{i, j)\ vs. the corre- 
sponding first 2000 coefficients of j)| with the largest 
modulus. S/N = 20 dB. 



Coefficient | DFT[ A(i,j) ] 

Fig. 9. Comparison of the results obtained by Tikhonov 
coupled with GCV in case of Poissonian and Gaussian 
noises, This figures correspond to the experiment shown 
in Fig. El Panel a) the coefficients \bg{i,j)\ and \bp{i,j)\ 
in decreasing order. The two horizontal lines represent the 
values of A for the two noises; b) corresponding GCV func- 
tions; c) coefficients and \bg{i,j)\ corresponding 
to the first 2000 coefficients of |y4(z,j)| shown in panel a). 
The vertical lines show the indices of closest to 
A. d) Ab{i,j) = \bg{i,j) - bp{i, j)\/\bg{i, j)\ vs. the corre- 
sponding first 2000 coefficients of j)| with the largest 
modulus. S/N = 30 dB. 
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Fig. 10. Comparison of the results obtained by Tikhonov 
coupled with GCV in case of Poissonian and Gaussian 
noises, This figures correspond to the experiment shown 
in Fig. 21 Panel a) the coefficients \bg{i,j)\ and \bp{i,j)\ 
in decreasing order. The two horizontal lines represent the 
values of A for the two noises; b) corresponding GCV func- 
tions; c) coefficients |6p(i,j)| and \bg(i,j)\ corresponding 
to the first 2000 coefficients of j)| shown in panel a). 
The vertical lines show the indices of |&(i,j)| closest to 
A. d) Ab{i,j) = \bg{i,j) - bp{i, j)\/\bg{i, j)\ vs. the corre- 
sponding first 2000 coefficients of with the largest 
modulus. S/N = 40 dB. 



Fig. 11. Comparison of the results obtained by Tikhonov 
coupled with GCV in case of Poissonian and Gaussian 
noises, This figures correspond to the experiment shown 
in Fig. El Panel a) the coefficients \bg{i,j)\ and \bp{i,j)\ 
in decreasing order. The two horizontal lines represent the 
values of A for the two noises; b) corresponding GCV func- 
tions; c) coefficients \bp{i,j)\ and \bg{i,j)\ corresponding 
to the first 2000 coefficients of j)| shown in panel a). 
The vertical lines show the indices of |6(i,_;/)| closest to 
A. d) Ab{i,j) = \bg{i,j)-bp{i,j)\/\bg{i,j)\ vs. the corre- 
sponding first 2000 coefficients of j)| with the largest 
modulus. S/N = 60 dB. 
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Fig. 12. LW deblurring of the images shown in Fig. |21 
rmsr and rms/ are the rms of the true residuals calculated 
on the entire image and only on the pixels corresponding 
to the satellite, respectively. 



Fig. 13. LW deblurring of the images shown in Fig. 13 
rmsT and rms/ are the rms of the true residuals calculated 
on the entire image and only on the pixels corresponding 
to the satellite, respectively. 
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Fig. 14. LW deblurring of the images shown in Fig. 0] Fig. 15. LW deblurring of the images shown in Fig. 

rmsT and rms/ are the rms of the true residuals calculated rmsr and rms/ are the rms of the true residuals calculated 

on the entire image and only on the pixels corresponding on the entire image and only on the pixels corresponding 

to the satellite, respectively. to the satellite, respectively. 
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Fig. 16. Convergence rate of the LW algorithm applied to pig. 18. Convergence rate of the LW algorithm apphed to 
the problem of deblurring the images shown in Fig. El the problem of deblurring the images shown in Fig. H 



Landweber Method - S/N = 30 dB 



Landweber Method - S/N = 60 dB 
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Fig. 17. Convergence rate of the LW algorithm applied to ^ig. 19. Convergence rate of the LW algorithm apphed to 
the problem of deblurring the images shown in Fig. ^he problem of deblurring the images shown in Fig. 
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Fig. 20. \\xk — a;||/||a;|| vs. the number of iterations. The 
object of interest, located in the center of the image, is 
a two-dimensional Gaussian with circular symmetry and 
dispersion (in pixels) given in the top of each panel. It is 
superimposed to a background whose level is 1% the peak 
value of the blurred image. The PSF is a two-dimensional 
Gaussian with a dispersion of 6 pixels. The size of the 
image are 256 x 256 pixels. The noise is Poissonian with 
peak S/N = 30 dB. Two deblurring algorithms are used: 
Richardson-Lucy (RL) and LW. 



Fig. 22. As in Fig. I^Olbut with S/N = 40 dB. 





Number of Iterations 



Number of Iterations 



Fig. 23. As in Fig. |2Dlbut with S/N = 40 dB. 
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Fig. 21. As in Fig.EHl 
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Fig. 24. ||a;fc — a;||/||a;|| vs. the number of iterations. The 
object of interest, located in the center of the image, is of a 
two-dimensional rectangular function, with side length (in 
pixels) given in the top of each panel. It is superimposed 
to a background whose level is 1% the peak value of the 
blurred image. The PSF is a two-dimensional Gaussian 
with a dispersion of 6 pixels. The size of the image are 
128 X 128 pixels. The noise is Poissonian with peak S/N = 
30 dB. Two deblurring algorithms are used: Richardson- 
Lucy (RL) and LW. 
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Fig. 26. As in Fig. Elbut with S/N = 40 dB. 
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Fig. 27. As in Fig.Elbut with S/N = 40 dB. 
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Fig. 25. As in Fig. El 
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Fig. 28. fjfc vs. the number of iterations for RL in the 
case of a Gaussian source with various values of ctq and a 
Gaussian PSF with dispersion (jp = 6. 
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Fig. 29. i?^/^((T^) vs. the number of iterations for the 
cases shown in Fig. ESI 
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